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ABSTRACT 


A  new  approach  to  the  application  of  maximum  entropy  principles  in  image  restoration  is 
being  developed.  It  is  an  iterative  process  which  combines  a  Wiener  filter  restoration 
and  modifications  to  the  wiener  image  that  are  Entropy-based.  The  result  is  an  image  with 
well-restored  frequency  content  and  very  little  of  the  spuriousness  commonly  introduced  by 
inverse  filters.  The  algorithm  is  fast,  stable  to  convergence,  and  will  accommodate  any 
specifiable  distorting  function. 


_ INTRODUCTION 


It  is  frequently  desirable  to  remove  taking-system  degradation  from  images  in  an 
attempt  to  obtain  a  truer  object  record.  The  most  conceptually  straightforward  of  these 
restoration  processes  are  linear  inversion  techniques,  in  which  the  distorting  spread 
function  is  effectively  deconvolved.  Often,  in  practice,  this  amounts  to  a  fre^ency 
domain  division.  While  the  simplicity  of  inversion  methods  is  attractive,  inherent  in  the 
frequency  domain  division  is  the  introduction  of  large  spurious  artifacts  in  the  regions 
where  the  MTF (divisor)  is  near  zero.  The  resulting  restorations  are  marked  by  large 
periodic  'side  lobes'  and  noise  that  are  soley  due  to  the  restoration  process,  and  do  not 
represent  the  object.  Such  results  are  usually  undesirable. 

Another  approach  to  object  estimation  from  images  is  probablistic  and  involves  the 
concept  of  entropy  in  the  information  capacity  sense  of  the  word.  Maximum  entropy 
(Maxent)  restoration  methods  use  the  gathered  image  and  any  other  available  information 
about  the  object  to  predict,  without  making  assumptions,  that  object  which  most  likely 
gave  rise  to  the  recorded  image.  Since  these  estimations  begin  with  no  assumptions  about 
the  object,  the  first  iteration  is  a  flat,  uniform-gray  field.  As  the  restoration 
proceeds,  a  likely  object  rises  out  of  the  flat  field  consistent  with  the  observed  data 
and  any  other  known  information.  After  this  process,  the  resulting  object  estimate 
consists  only  of  information  that  was  present  in  the  image  data.  No  spurious  noise  or 
artifacts  are  introduced  by  the  process. 

We  present  an  image  restoration  algorithm  which  is  a  hybrid  of  a  linear  inversion 
technique  and  the  Maxent  concept.  The  results  reported  here  are  an  extension  of  research 
reported  in  references  1  and  2. 


_ PEBIYATIOW 

Assume  a  subject  image  which  is  digital  and  was  formed  by  the  process 

y(k)  -  x(k)  *  p(k)  t-  n(k) , 

where  denotes  convolution,  y  is  the  observed  data,  x  is  the  object,  p  is  the  imaging 
system  point  spread  function,  and  n  is  noise.  At  least  y  and  p  are  known. 

The  problem  is  to  estimate  x(k)  from  restoration  of  y(k)  .  Once  an  estimate  of  x  is 
arrived  at,  x,  is  convolved  with  p  to  form  an  estimated  noiseless  image,  g,(h): 
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ggCk)  »  Xg(k)  *  P(k) 

From  this  quantity,  the  mean  square  error  is  formed: 

MSE  ■  I  [y(k) -ge(k) 

The  entropy  of  Xe(k)  is  defined  as 


S  -  -r  Xe(k)ln[Xg(k) ]. 

In  many  Maxent  approaches^ ' ^ ®  entropy  is  maximized  without  letting  g^Ck)  depart 
significantly  from  the  observed  data,  y(k) .  In  particular,  a  quantity,  Q,  is  formed 
which,  by  use  of  a  Lagrange  multiplier,  binds  entropy  and  MSE  together  such  that  a  trade¬ 
off  is  made  between  them  when  an  object  estimate  is  formed.  To  bring  the  solution  to  an 
differentiate  Q,  set  the  result  to  zero,  and  solve  for  x_(k).  The  resulting 
implicit  equation 


Xg(k)  »  exp{-l  +  Xty(k)  -  g^Ck)]  *  p(-k)), 
is  the  basis  for  an  iterative  solution, 

Xe"‘^^(k)  -  exp(-l  X[y(k)  -  g^^Ck)  ]  *  p(-k)), 
in  which  XgQ(k)  is  a  constant. 

Restorations  from  this  kind  of  iterative  process  are  smooth,  noise-free,  and  generally 
result  in  a  2-4x  resolution  increase  after  many  tens  of  iterations.  It  should  be  noted 
that  in  this  form,  the  solution  tends  to  be  unstable,  typically  requiring  the  "slow" 
averaging  in  of  new  iterations;  e.g. 

Xg"(k)  *  0.95(Xg""^{k))  +  0.05(Xcg"(k)) , 

where  Xgg"(k)  is  the  estimate  that  would  be  directly  calculated  by  the  algorithm.  It  is 
this  characteristic  that  both  allows  for  eventual  convergence  and  leads  to  the  possible 
need  for  a  hundred  or  more  iterations. 

It  is  desirable  to  find  a  solution  which  results  in  a  greater  resolution  increase, 
converges  more  quickly,  is  more  stable,  and  less  compute-intensive,  especially  when  one  is 
faced  with  the  restoration  of  a  large  database.  It  is  well  known  that  inversion 
techniques  do  an  adequate  job  of  resolution  restoration,  and  that  the  maximization  of 
entropy  in  a  solution  formulation  tends  toward  smooth,  flat  object  estimates.  It  is  a 
reasonable  idea  then,  to  use  an  inverse  filter  to  do  the  bulk  of  the  resolution 
restoration,  and  apply  entropy  constraints  in  an  attempt  to  defeat  the  inversion- induced 
artifacts. 

If,  in  the  above  derivation,  after  the  differentiation  of  Q,  one  solves  for  x_  by 
taking  the  Fourier  Transform  of  both  sides  instead  of  exponentiating,  a  solution  with  two 
terms  results.  The  first  is  a  wiener  filter,  and  the  second  is  an  entropy  gradient. 

Here  entropy  is  defined  as  S  -  -J  h(k) InCh(k) ] ,  where  h(k)  -  Xg(k)/E  Xg(k). 

Differentiate  Q  with  respect  to  x(k)  and  set  the  result  to  zero: 

Q  -  S  -  v/2  MSE 

-S  -  log  h(k)  +  V  [y(k)  -  g(k)]  *  p(-k)  -  0 
Fourier  Transform  both  sides  and  solve  for  X  in  the  frequency  domain: 

X(f)  -  W(f)  -  V  B(f) 

W(f)  is  a  Wiener  filter: 


and  B(f)  is 


W(f)  -  y(f)  p*(f) 

■|P(f) +  C 

B(f)  -  D(f) 

|P(f) I*  ♦  C 
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where 


d(lc)  =  S  +  log  h(k)  . 


Note  that  b(lc)  is  a  filtered  entropy  gradient  since  the  entropy  term,  d(k) ,  is  filtered  by 
1/!P|^  to  fora  b'(k)  . 


A  small  constant,  C,  is  added  to  the  denominator  to  control  the  result  for  very  small 
values  of  P(f).  This  is  usual  in  Wiener  inversions.  C  is  an  expression  of  noise:signal 
and  is  treated  as  a  constant,  though  it  could  be  functional. 

Choice  of  the  initial  value  of  C  and  the  intermediate  values  of  the  entropy  gradient 
weighting,  v,  are  important  for  stability  and  convergence.  As  C  goes  to  zero,  the 
starting  estimate  for  the  entropy  iterations  becomes  sharper,  so  one  wants  to  choose  C 
such  that  the  Wiener  filter  performs  the  bulk  of  the  possible  resolution  increase,  but 
without  making  the  starting  estimate  so  oscillatory  that  the  algorithm  will  never  converge 
to  a  solution.  Currently  an  attempt  is  made  to  choose  Cq  to  minimize  noise: signal  in  the 
Wiener  inversion.  This  method  of  defining  was  chosen  because  it  is  based  only  on  the 
data  (no  user  experience  required)  ,  it  allows  C  to  be  some  real  representation  of  its 
theoretical  function,  and  it  can  be  done  automatically  by  the  algorithm.  This  method 
works  well  for  most,  but  not  all,  of  the  images  we  have  restored,  so  additional  ways  to 
choose  Cq  are  being  investigated. 

All  of  the  finesse  of  this  restoration  algorithm  lies  in  the  iterative  application  of 
the  entropy  gradient.  The  shape  of  the  gradient  is  given  unambiguously  as  the  second  term 
of  the  solution  for  x(k)  by  the  derivation.  Its  weighting,  however,  remains  to  be  defined 
at  each  iteration.  Since  the  algorithm  is  a  spectral  subtraction  of  a  weighted  entropy 
gradient  from  a  Wiener  image,  the  weighting  must  be  large  enough  to  be  effective  against 
the  Wiener— induced  artifacts,  but  not  so  large  as  to  introduce  additional  artifacts. 
Several  calculable  attributos  of  each  iteration,  such  as  entropy,  MSE,  noise,  and  contrast 
are  used  at  each  iteration  to  choose  a  weight  which  will  drive  the  object  estimate  closer 
to  a  convergent  solution. 


3 .  RESULTS 

A  one-dimensional  example  demonstra -es  this  restoration  concept  graphically.  Figure  l 
is  a  cut  through  a  computer-simulated  object.  Figure  2  is  the  sane  image  degraded  by  21- 
pixel  linear  motion.  Figure  3  is  a  graphical  representation  of  the  two  terms  which 
comprise  the  restoration;  3a  is  the  Wiener  inversion  (notice  the  periodic  side  lobes  that 
have  been  introduced) ,  and  3b  is  the  entropy  gradient  term  for  the  first  iteration. 
Clearly,  the  subtraction  of  3b  from  3a  mainly  affects  (reduces)  the  spurious  side  lobes. 
The  effect  of  the  spectral  subtraction  is  shown  in  Fi^re  4a.  The  original  object  is 
obviously  restored,  but  still  somewhat  degraded  with  noise.  Figure  4b  shows  the  entropy 
gradient  for  the  next  iteration.  Notice  that  it  is  smaller  in  amplitude,  in  accordance 
with  the  reduced  noise  in  the  estimate  it  must  currently  act  upon. 


Restoration  of  a  two-dimensional  simulated  point  source  group  tests  how  well  the 
algorithm  preserves  relative  intensities,  as  well  as  showing  up  _  any  characteristic 
deformations  the  algorithm  might  introduce.  The  object  is  the  four  point  sources  shown  in 
Figure  5.  In  Figure  6,  that  object  has  been  degraded  by  convolution  with  a  real  point 
spread  function  from  the  IRAS  (Infrared  Astronomical  Satellite)  telescope.  The  detectors 
in  the  IRAS  focal  plane  are  rectangular,  have  asymmetric  responses  in  both  dimensions,  and 
are  very  large  compared  to  most  of  the  celestial  objects  imaged  onto  them.  The  IRAS  psf, 
then,  is  mainly  the  detector  response  profile.  The  simulated  IRAS  image  in  Figure  6  has 
had  a  of  peak,  Gaussian  noise  added  to  it.  Figure  7  shows  the  restoration  of  the  point 
source  group  after  20  iterations  of  the  FIER  (Filtered  Entropy  Restoration)  algorithm. 
Both  the  relative  and  absolute  peak  intensities  have  been  quite  well  restored.  The 
residual  noise  is  on  the  order  of  2%. 


Among  the  IRAS-scanned  objects  is  a  double  galaxy,  N2992,  which  is  effectively  a  pair 
of  point  sources  relative  to  IRAS  sampling  rates.  Figure  8  shows  data  from  the  60  urn  band 
IRAS  coverage  of  N2992.  The  system  psf  shape  is  obvious  in  this  data.  In  Figure  9,  the 
Wiener  inversion  is  shown.  The  two  point  sources  are  immediately  resolved  by  this  first 
step  of  FIER.  However,  many  artifacts  have  been  introduced.  After  approximately  20 
iterations  of  application  of  the  entropy  gradient,  a  restored  image,  or  object  estimate 
results  which  contains  the  full  wiener  resolution  restoration  in  the  object  region,  and  is 
nearly  free  of  side  lobes  and  other  spurious  noise.  The  small  side  lobes  which  do  remain 
are  on  the  order  of  1%  of  peak  intensity.  See  Figure  10. 


□ 
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Avoi  itibiiity  Codes 


An  image  restoration  algorithm  which  combines  the  most  desirable  properties  of  both 
linear  inversion  techniques  and  Entropy-based  estimation  techniques  has  been  developed  and 
used  to  restore  astronomical  images.  This  filtered  entropy  algorithm  uses  a  Wiener  filter 
to  effect  resolution  restoration.  It  then  iteratively  applies  an  entropy-based  gradient 
in  the  frequency  domain  to  defeat  the  spuriousness  generated  by  the  inversion.  The 
algorithm  is  stable,  converges  quickly  without  averaging  with  previous  iterations,  and 
requires  no  assumptions  about  either  the  noise  in  the  observed  data,  or  the  shape  of  the 
psf;  any  specifiable  psf  can  be  used.  The  resulting  restorations,  or  object  estimations, 
for  objects  rising  out  of  a  low  intensity  surround  represent  a  6-7x  resolution  increase 
over  the  observed  data,  and  the  process-induced  residuals  are  approximately  1%  of  peak 
intensity  of  the  final  restored  image. 
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Figure  2.  One  dimensional  object  degraded  by  21-pixel  linear  motion 
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Figure  3.  Inverse  (Wiener)  filter  solution  and  entropy  gradient  term 
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Figure  7.  simulated  point  source  image  restored  by  FIER 
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Figur*  10.  IRAS  «0um  N2992  r«*tor«d  by  flER 


